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Heat  and  Water  Vapor  Flux  in  Conifer  Forest  Litter  and  Duff: 
A  Theoretical  Model 


Michael  A.  Fosberg 


Introduction 

The  uppermost  layers  of  the  forest  floor  are  a  com- 
plex interactive  component  of  the  natural  water  bal- 
ance. The  organic  particles  found  in  these  layers  have 
water  sorption  properties  that  depend  on  their  previous 
state  and  on  the  chemistry  and  cellular  structure  of  the 
parent  organic  material.  The  physical  structure  of  the 
layers  depends  also  on  the  microclimate  of  the  forest 
stands,  as  well  as  the  physiology-dependent  organic 
decay  rates  of  the  particles. 

These  particles  are  rarely  in  equilibrium  with  the 
microenvironment  of  the  soil  voids,  and  almost  never 
in  equilibrium  with  the  external  environmental 
stresses  imposed  by  the  atmosphere.  Organic  soil 
particles  respond  to  the  soil  void  microenvironment 
more  slowly  than  changes  in  that  void  environment. 
An  individual  organic  soil  particle  is  both  a  source  and 
sink  for  water  and  energy  in  a  complex,  dynamic 
system.  Latent  and  sensible  heat  exchange  and  mass 
conservation  of  water  take  place  simultaneously  and 
interact  between  the  particles,  soil  voids,  and  the  ex- 
ternal environment.  These  processes  must  be  ac- 
counted for  through  an  interactive  model  to  physically 
depict  the  processes  in  accurate  detail. 

Wildfire  behavior  and  damage  are  strongly  in- 
fluenced by  the  moisture  content  of  the  organic  forest 
soils.  Prescribed  fire  for  removal  of  the  litter  and  duff 
for  seedbed  preparation,  and  for  hazard  reduction 
through  removal  of  the  understory  also  depends  on  the 
litter  and  duff  moisture  content.  Reliable  means  of 
predicting  the  moisture  contents  of  litter  and  duff  are 
required  for  predicting  effort  required  in  wildfire  con- 
trol and  for  safe,  beneficial  use  of  prescribed  fire. 

In  this  Paper,  a  theoretical  model  based  on  conser- 
vation of  energy  and  mass  was  developed  for  predicting 
the  moisture  content  of  litter  and  duff  under  high 
moisture  stress.  The  prediction  model  was  obtained  by 
first  considering  restrictive  analytical  solutions  of 
simple  functions  for  diffusion  of  single  layers  of  porous 
media.  These  simple  solutions  for  fixed  boundary  con- 
ditions defined  a  general  equation  to  evaluate  timelags 
from  numerical  solutions  of  nonlinear  differential 
equations.  The  bovmdary  conditions  were  then  relaxed 
and  allowed  to  vary  as  small  time-dependent  step  func- 
tions in  a  second  analytical  solution.  This  solution 
defined  a  more  general  case,  and  provided  the  bases  for 


numerical  experiments  with  the  nonlinear  equations 
with  variable  boundary  conditions  to  evaluate  and 
calibrate  the  coefficients  of  the  solution.  This  general 
solution  represented  the  response  of  the  litter  and  duff 
to  natural  environmental  stresses,  and  provided  the 
functional  framework  for  predicting  the  moisture  and 
temperature  profiles. 

Analysis  of  Previous  Work 

Soil  moisture  and  heat  transfer  analysis  has 
traditionally  been  applied  to  saturated  or  partially 
saturated  soils.  This  emphasis  on  the  presence  of  liquid- 
water  transfer  stems  from  agricultural  needs  for  mois- 
ture in  the  root  zone,  and  for  conservation  of  impounded 
water.  While  the  complete  processes  of  heat  and  mass 
transfer  between  soil  particles  and  soil  voids  and 
between  adjacent  voids  are  qualitatively  known,  most 
previous  quantitative  work  has  concentrated  on  the 
processes  associated  with  liquid  water  (Chudnovskii 
1962,  Childs  1969).  Under  these  conditions,  the  largest 
mass  of  soil  water  is  in  the  soil  voids,  and  the  moisture 
bound  within  and  on  the  soil  particles  contributes  little 
to  the  water  balance  as  sources  or  sinks  of  moisture.  Be- 
cause of  this  emphasis  on  liquid  water,  most  analytical 
and  theoretical  work  has  been  based  on  a  simple  Darcy 
type  of  flow  where  the  soil  structure  can  be  accounted  for 
by  simple  diffusivity  or  hydraulic  conductivity.  The 
assumption  that  the  particles  do  not  contribute  mater- 
ially to  the  water  and  energy  budgets  breaks  down  if  no 
liquid  water  is  present,  because  the  total  heat  content 
and  moisture  in  the  voids  is  small  compared  to  those  in 
the  particles  (USDA-FPL  1974,  Jackson  1964c). 

Studies  of  soil  aeration  with  nonsorbing  porous 
media  demonstrated  that  the  diffusion  rates  of  gases 
and  vapors  are  a  substantial  fraction  of  the  molecular 
diffusion  rates  in  free  air.  The  diffusion  of  gases  in 
porous  media  may  be  as  high  as  10  percent  of  that  in 
free  air,  even  with  porosities  as  low  as  15  percent 
(Penman  1940;  Taylor  1949;  Millington  1959;  Currie 
1960a,  19606,  1961;  Millington  and  Schearer  1971; 
Ayreset  al.  1972).  For  porosities  of  70  to  90  percent,  the 
range  common  to  duff  and  litter,  the  diffusion  rate  is  60 
to  80  percent  of  that  in  the  free  air.  Thus,  soils  in 
general  and  litter  and  duff  horizons  in  particular  do  not 
markedly  reduce  the  gas  exchange  by  their  structure. 
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Experimental  determinations  of  soil-void  vapor 
diffusivities  with  sorbing  materials  have  produced 
mixed  results;  studies  by  de  Jong  and  Schappert  (1972) 
and  Hanks  (1958)  show  agreement  with  the  nonsorbing 
materials.  Jackson  (1964a,  19646)  calculated  diffu- 
sivities two  to  three  orders  of  magnitude  lower  than 
those  found  in  nonsorbing  experiments.  Jackson's 
(1964c)  steady  state  experiments  provide  a  clue  on  how 
to  reconcile  this  discrepancy.  The  sorption  isotherms 
indicate  that  the  soil  particles  contain  substantial 
amounts  of  water  in  comparison  with  that  contained 
in  the  voids  at  equilibrium.  Jackson  assumed  that  the 
particles  are  always  in  equilibrium  with  the  voids.  If 
the  exchange  of  moisture  between  the  voids  and  part- 
icles is  time  dependent  (lags  behind),  then  use  of  the 
equilibrium  isotherms  to  evaluate  the  diffusivity  will 
underestiniate  the  diffusivity  in  the  voids.  If  the  part- 
icles respond  very  slowly  to  stress,  then  the  estimate  of 
diffusivity  will  be  very  low.  Only  as  their  response  rate 
approaches  the  instantaneous  will  the  equilibrium 
isotherm  technique  be  valid.  Therefore,  one  must 
know  the  response  rates  of  individual  particles  as  well 
as  the  bulk  diffusivities  of  the  soil  matrix  to  calculate 
the  heat  and  mass  exchange. 

Two  sets  of  characteristics  of  particle  behavior 
need  to  be  defined  to  properly  account  for  heat  and  mass 
transfer  between  the  particles  and  voids:  timelags  and 
equilibrium  isotherms.  The  timelag  is  defined  through 
the  time  interval  during  which  a  particle  will  gain  or 
lose  heat  or  moisture  from  a  homogeneous  state  after 
a  specific  change  in  the  environment.  The  change  takes 
place  as  an  exponential  decay  to  the  new  equilibrium 
condition.  The  timelag  is  then  defined  as  the  time  inter- 
val required  to  achieve  1-1/e  (approximately  62  per- 
cent) of  the  potential  change.  Theoretical  analysis 
under  laboratory  conditions  has  defined  the  timelag  in 
terms  of  particle  properties  (Byram,^  Fosberg  1970). 
This  theoretical  work  was  extended  to  field  situations 
by  generalizing  the  boundary  conditions  in  the  solu- 
tions (Fosberg  1972,  1973).  The  equations  in  Fosberg 
(1972)  provide  the  basis  for  moisture  exchange  in  the 
soil  particles,  and  the  equations  in  Fosberg  (1973)  do 
the  same  for  heat  exchange.  Moisture  timelags  for 
many  of  the  litter  and  duff  particles  have  been  analyzed 
by  Nelson  (1969),  Van  Wagner  (1969),  Fosberg  (1970), 
and  Mutch  and  Gastineau  (1970).  Data  on  thermal 
timelags  are  far  less  available.  The  only  compilation 
for  organic  materials  is  in  Fosberg  (1973).  Unfortu- 
nately, no  data  exist  for  either  moisture  or  thermal 
timelags  of  soil  particles. 

The  second  set  of  characteristics  that  must  be 
known  are  the  equilibrium  isotherms  for  moisture  and 
specific  heat  of  the  soil,  litter,  and  duff  particles.  These 
isotherms  define  the  boundary  conditions  required  for 
heat  and  water  exchange.  The  equilibrium  moisture  iso- 

^Byram,  George  M.  1963.  An  analysis  of  the  drying  proc- 
ess in  forest  fuel  materials.  Paper  presented  at  the  1963 
International  Symposium  on  humidity  and  moisture.  Wash. 
D.C.  May  20-23,  1963.  38  p. 


therms  for  soil  particles  have  been  compiled  by  Baver 
(1956),  Jackson  (1964c),  and  Miyamoto  et  at.  (1972). 
As  with  the  timelags,  there  are  numerous  sources  of 
equilibrium  isotherm  data  for  woody  materials.  The 
USD  A  Forest  Products  Laboratory  (1974)  has  compiled 
a  comprehensive  set  of  tables  for  commercial  woods. 
Data  for  fine  particles  found  in  litter  and  duff  are  pre- 
sented in  Mutch  and  Gastineau  (1970)  and  Blackmarr 
(1971).  The  thermal  properties  of  soil  particles  are  well 
documented  in  soil  texts  (Baver  1956,  Chudnovaskii 
1962)  and  those  of  woody  materials  in  Siau  (1971). 

Analysis  of  the  soil  physics  and  fuel  moisture  liter- 
ature lead  directly  to  a  conceptual  model  of  energy  and 
vapor  flow  in  forest  litter  and  duff.  The  particles  respond 
slowly  to  changes  in  the  void  microenvironment,  while 
the  environment  of  the  voids  can  change  rapidly.  This 
process  was  recognized  near  the  beginning  of  fuel 
moisture  research.  Early  papers  clearly  defined  the 
litter  and  duff  as  well  aerated  (Gast  and  Stickel  1929) 
and  stated  that  the  particle  response  lagged  behind 
changes  in  humidity  (Harkness  1939).  This  slow  re- 
sponse of  litter  and  duff  particles  to  liquid  water  was 
observed  more  recently  by  Clary  and  Ffolliott  (1969). 
Most  of  the  recent  work  on  duff  and  litter  moisture, 
particularly  that  intended  for  predictive  methods,  has 
departed  from  relations  depicting  the  physical  proc- 
esses, and  have  concentrated  on  developing  a  drought 
index  or  build-up  index.  The  approach  with  the  syn- 
thetic indices  required  regression  equations  to  relate 
the  indices  to  duff  and  litter  (Jarvis  and  Tucker  1968, 
Johnson  1968,  Muraro  and  Lawson  1970).  While  this 
approach  provides  easily  used  tools,  the  indices  are  not 
directly  related  to  duff  and  litter  moisture,  and  a  new 
set  of  regressions  must  be  obtained  for  each  site.  Van 
Wagner  (1970)  approached  the  prediction  problem  by  a 
method  more  closely  related  to  the  physical  processes. 
Because  he  used  a  mathematical  model  developed  for 
fixed  boundary  conditions,  his  model  does  not  fully 
depict  the  physical  processes.  Van  Wagner  transformed 
the  model  from  a  constant  environment  solution  to  a 
field  solution  through  a  variable  timelag.  This  variable 
timelag  combined  both  the  lag  produced  by  the  physical 
structure  and  the  variations  in  moisture  change  pro- 
duced by  a  varying  environment. 


Structure  of  the  Model 

Four  distinct  processes  take  place  simultaneously 
in  the  litter  and  duff  system.  The  organic  particles  ex- 
change heat  and  water  vapor  with  the  surrounding 
voids.  The  exchange  rates  are  determined  by  (a)  the 
boundary  conditions  imposed  by  the  temperature  and 
relative  humidity  in  the  voids,  (b)  the  excess  or  deficit 
of  heat  and  moisture  in  the  particles  in  comparison  to 
the  boundary  condition  imposed  by  the  void  atmo- 
sphere, and  (c)  the  rate  at  which  the  particles  can  take 
up  or  give  off  heat  and  water  vapor.  The  last  two  pro- 
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cesses  involve  exchange  of  heat  and  water  vapor  be- 
tween adjacent  voids  in  response  to  the  concentration 
gradients. 

The  processes  at  a  given  level  in  the  system  are 
illustrated  schematically  in  figure  1.  The  four  corners 
represent  the  state  variables  of  heat  and  vapor  in  the 
particles  and  in  the  voids.  The  bold,  horizontal  double- 
headed  arrovvs  define  the  dynamic  processes  of  ex- 
change between  the  particles  and  the  voids.  The  bold 
vertical  arrows  indicate  the  exchange  between  ad- 
jacent voids.  The  light  arrows  and  interior  cells  define 
the  flow  of  energy  and  vapor  and  the  interactive  process 
of  exchange  between  the  particles  and  voids.  The  ab- 
solute moisture  content  and  temperature  of  ithe  voids 
specify  the  relative  humidity  of  the  voids.  The  void 
relative  humidity  in  turn  specifies  the  equilibrium 
moisture  content  that  the  particles  would  eventually 
achieve.  This  equilibrium  moisture  content  represents, 
then,  the  environment  to  which  the  particles  must 
respond,  and  is  expressed  as  the  boundary  condition 
from  the  vapor  transfer  in  to  or  out  of  the  particles. 

The  transfer  of  vapor  in  to  or  out  of  a  particle  de- 
pends on  the  magnitude  and  direction  of  the  difference 
between  the  actual  particle  moisture  content  and  its 
equilibrium  moisture  content.  The  particles  do  not 
respond  instantaneously  to  changes  in  the  environ- 
ment. The  rate  of  response  is  expressed  through  the 
timelag  in  a  quasi-exponential  approach  to  the  equi- 


librium. The  particle  timelag  contains  the  information 
of  size,  internal  moisture  diffusivity,  and  shape  of  the 
particle.  Heat  transfer  proceeds  in  a  similar  manner. 

The  temperature  and  heat  content  of  the  voids 
define  the  boundary  condition  imposed  on  the  particles. 
Like  the  moisture  exchange  between  the  particles  and 
voids,  heat  transfer  depends  on  the  temperature  grad- 
ient and  on  the  physical  characteristics  of  size,  shape, 
and  specific  heat.  In  addition,  constraints  are  imposed 
by  the  influence  of  moisture  on  the  thermal  properties 
and  by  the  second  law  of  thermodynamics  because  of 
the  large  difference  in  specific  heats.  The  particle 
moisture  content  influences  the  heat  exchange  through 
effects  on  the  specific  heat  and  through  volumetric 
changes  of  the  particles.  Also,  because  of  the  large 
differences  in  specific  heats  and  because  heat  cannot 
be  moved  against  a  temperature  gradient,  the  laws  of 
conservation  of  heat  and  mass  governing  the  processes 
must  be  supplemented  with  the  second  law  of  thermo- 
dynamics. 

The  Mathematical  Model 

This  descriptive  model  of  heat  and  vapor  transfer 
in  duff  and  litter  is  based  on  two  underlying  concepts: 
(1)  diffusive  exchange  between  voids  is  rapid,  and  (2) 
exchange  between  voids  and  particles  is  slow,  and  is  the 
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Figure  1  .—Schematic  diagram  of  heat  and  mass  transfer  processes  in  the  internal  layers  of  conifer  forest  litter 

and  duff. 
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limiting  factor  in  heat  and  vapor  exchange.  A  theoret- 
ically derived  model  can  now  be  constructed  from  these 
descriptive  concepts. 

The  transport  of  v/ater  vapor  through  the  voids  is 
governed  by  the  conservation  of  mass  equations,  where 
changes  in  density  depend  on  diffusion  of  mass  and  on 
the  vapor  taken  up  or  given  off  by  the  individual  part- 
icles. The  continuity  equation  for  vertical  transport 
in  a  horizontally  homogeneous  layer  is: 


3p 
 i 

at 


3^P        .  3p 
v  ,    8v  V 


3z' 


3z  3z 


±  X 


appropriate  provided  t  is  replaced  by  8t,  the  finite  time 
step.  The  continuity  equation  for  water  vapor  is, 
therefore: 


9p  9p 
v  ^  Y  

9t  9z  9z 


9v  v 


-  X 
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where  x  is  approximated  over  a  small  finite  interval  by 
6mp_ 


X  = 


B 


[2] 


where  p^is  the  density  of  water  vapor  in  the  voids,  v  is 
the  molecular  diffusivity  governing  the  transport  be- 
tween voids,  and  x  is  the  rate  of  exchange  of  water  vap- 
or between  the  particles  and  the  voids.  The  moisture 
content  of  the  litter  and  duff  particles  is  expressed  on 
an  ovendry  gravimetric  base.  Gain  or  loss  of  moisture 
content  by  the  particles  must  be  transformed  to  units 
consistent  with  those  of  the  continuity  equation.  The 
moisture  content  of  the  particle  is  expressed  as  the 
mass  of  water  divided  by  the  ovendry  mass  of  the 
particle,  and  the  vapor  density  is  expressed  as  the  mass 
of  water  vapor  in  the  void  volume.  The  particle  moisture 
content  is  transformed  to  equivalent  vapor  density 
through  multiplication  by  the  bulk  density  of  the  layer 
and  division  by  the  porosity: 


where  m  is  the  particle  gravimetric  content,  Pb  is  the 
layer  bulk  density,  and  4>  is  the  porosity.  The  relation- 
ship between  particle  changes  in  moisture  and  the  cor- 
responding moisture  taken  up  or  given  off  to  the  voids 
is  now  obtained  by  differentiation  and  thus  defines  the 
source  or  sink  term  x- 


and 


dm 
dt 


dp' 


dt 


=  X 


The  changes  in  particle  moisture  content  produced 
by  environmental  stress  derived  by  Fosberg  (1972)  are: 


6m      T  ^ 

-T—  =  1  -  C  exp 

Am  m 


t 
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where  5m  is  the  actual  moisture  exchange.  Am  is  the 
potential  exchange  if  the  particle  responded  instantan- 
eously, {  is  a  coefficient  reflecting  the  internal  distribu- 
tion of  moisture,  t  j^is  the  particle  timelag,  and  t  is 
the  time  over  which  the  environmental  stress  acts.  The 
particle  moisture  exchange  may  be  differentiated  and 
substituted  directly  into  the  continuity  equation.  Since 
the  continuity  equation  is  nonlinear,  it  will  be  solved 
by  finite  difference  methods;  this  form  of  the  equation  is 
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[3] 


The  transfer  of  heat  through  the  voids  is  governed 
by  the  first  law  of  thermodynamics.  Since  the  partic- 
ular problem  is  concerned  with  a  shallow  system 
vertically,  the  first  law  may  be  approximated  as 


dT 
'p  dt 


dt 


dT 
dt 


where  Cp  is  the  specific  heat  at  constant  pressure,  T  is 
the  temperature,  and  s  is  the  source  or  sink  of  heat  pro- 
vided by  the  particles.  As  with  the  water  vapor,  the 
heat  is  exchanged  between  voids  by  diffusion.  The  first 
law  of  thermodynamics  is: 


:  9T 
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where  k  is  the  molecular  diffusivity  in  the  voids. 

The  thermal  source  or  sink  of  heat  provided  by  the 
litter  and  duff  particles  is  also  obtained  from  the  first 
law  of  thermodynamics.  Since  energy  must  be  con- 
served, any  exchange  of  heat  must  be  equal. 


6H. 


6H 


The  total  heat  of  the  air  is: 

6H    =  M    c  6T 
a       a    p  a 

where  is  the  heat  content  of  the  voids,  is  the 
mass  of  air  in  the  voids,  and  T^  is  the  void  air  temper- 
ature. Similarly,  for  the  particles: 


6h. 


and  the  void  heat  change  is: 
^F 


c  6t  = 
p  a 


M 


4 


r 


Now,  by  multiplying  this  expression  by  a  sequence 
of  unity  ratios: 


c  6T  = 
P  a 


mm 


^F  (1  - 


The  heat  exchange  is  expressed  in  terms  of  readily 
available  properties.  The  V's  refer  to  volume  and  the 
subscripts  F,  V,  and  B  represent  particle,  void,  and 
total  bed,  respectively. 

The  dynamic  response  of  temperature  changes  in 
the  particles  is  given  in  Fosberg  (1973)  and  is: 

6T  T  ^  -t 
^  =  1  -       exp  - 

where  ST  is  the  actual  temperature  change,  and  AT 
is  the  potential  change  for  an  instantaneous  response, 
Cj  is  the  internal  distribution  coefficient,  and  is  the 
thermal  response  time  or  timelag. 

Because  pp  cp  is  much  greater  than  Cp,  the  sec- 
ond law  of  thermodynamics  must  also  be  considered. 
It  is  not  necessary  to  quantitatively  consider  the  second 
law  in  the  numerical  solutions,  however.  If  the  time 
step  is  selected  sufficiently  small,  so  that  the  heat  flow 
from  the  high  temperature  to  the  lower  temperature 
is  maintained  in  the  finite  difference  calculation  ,  it 
is  sufficient  to  calculate  the  heat  exchange  based  only 
on  the  first  law.  The  equation  for  temperature  change 
in  the  voids  is: 


9t         9z2      3z  9z  c 


[4] 


The  temperature  change  in  the  voids  by  exchange  with 
the  particles  is  then 


=  6T. 


F  p  * 


[5] 


[3]  and  the  gain  or  loss  of  temperature  is  given  by 
equation  [5].  The  interface  of  mass  exchange  between 
equation  [1]  and  equation  [3]  is  given  by  equation 
[2].  The  conservation  of  heat  exchange  with  the 
restriction  of  preserving  the  second  law  of  thermo- 
dynamics is  given  in  equation  [5].  With  these  six 
equations,  the  moisture  and  temperature  profile  dyna- 
mics are  fully  defined. 

These  equations  must  be  supplemented  by  a  set 
of  diagnostic  state  equations  to  fully  close  the  system. 
These  diagnostic  equations  relate  the  void  air  thermo- 
dynamics properties  to  the  void  stresses  imposed  by 
relative  humidity  on  the  moisture  and  the  heat  con- 
tent of  the  particles.  The  boundary  condition  of  mois- 
ture transfer  is  defined  as  the  equilibrium  moisture 
content.  This  equilibrium  moisture  content  is  strongly 
dependent  on  relative  humidity.  Since  the  temperature 
and  vapor  density  are  known  at  any  given  time  and 
point,  the  equations  of  state 


p  R  T 
v  V 


[7] 


can  be  used  to  calculate  the  vapor  pressure  in  the  voids. 
At  constant  pressure,  the  saturation  vapor  pressure  is 
a  function  of  temperature  only.  The  saturation  vapor 
pressure  is  calculated  from  the  Goff-Gratch  form- 
ulation presented  in  the  Smithsonian  meteorological 
tables  (List  1966)  and  is  reproduced  here  for  com- 
pleteness. The  saturation  vapor  pressure  is 


+  5.02808  log 


1.3816  X  10-7    10^1. 344  fl  -^1 


[8] 


-t-  8.1328  X  10"^    do  3.4919  11     „  |  l] 


and  the  temperature  change  in  the  particles  is 


«Tj,.Tj,  -  Tj,.. 


where  Tg  is  the  boiling  point  temperature  (373.16°K) 
and  e^  is  the  saturation  pressure  of  ordinary  water  at 
the  boiling  point  temperature  (1.013246  x  10  dynes 
cm'^i.  The  relative  humidity  (h)  is  then  specified  by: 


Equations  [1]  through  [6]  dcfiine  the  four  dy- 
namic processes  of  heat  and  vapor  change  and  the  two 
sets  of  relations  between  void  and  particle  exchange. 
The  equation  of  continuity  [equation  1]  and  the  first 
law  of  thermodynamics  [equation  4]  define  the  void- 
to-void  exchange  in  terms  of  void  properties  and  gains 
from  or  losses  to  the  individual  particles.  The  gain  or 
loss  of  water  vapor  by  the  particles  is  given  by  equation 


h  = 


lOOe 


[9] 


The  boundary  condition  required  for  vapor  flux 
to  and  from  the  particles  is  expressed  by  the  equili- 
brium moisture  content.  Values  of  the  equilibrium 
moisture  content  of  softwoods  as  a  function  of  temper- 
ature and  relative  humidity  are  given  in  the  Wood 
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Handbook  (USDA-FPL  1974).  Equilibrium  isotherms 
for  pine  needles  were  obtained  from  unpublished  re- 
sults of  Blackmarr."* 

The  only  other  function  needed  to  close  the  set  of 
equations  is  the  relation  for  specific  heat  of  the 
particles.  The  specific  heat  according  to  Siau  (1971)  is: 


m  +  0.324 

1  +  m 


cal  gm 


-1  °C-1 


[10] 


The  diffusion  of  water  vapor  or  heat  through  the 
voids  is  predicted  from  Millington  and  Schearer's 
(1971)  tortuosity  factor.  They  determined  that  the 
vapor  and  thermal  diffusivities  would  be  reduced  from 
the  free  air  values  by  the  ratios: 


2x 


[11a] 


develop  the  similarity  relationship  came  from  computer 
simulations.  This  model  provided  four  equations.  The 
first  two  predicted  the  time-dependent  behavior  of 
heat  and  vapor  change  in  the  porous  bed,  and  the  second 
two  provided  the  heat  and  mass  timelags  from  physical 
structure. 

The  second  model  was  developed  to  transform  the 
laboratory-type  solution  to  field  situations  where  the 
boundary  conditions  are  continually  changing.  As  in 
the  laboratory  model,  this  model  was  developed  from 
analytical  and  numerical  solutions.  The  analytical 
solution  was  obtained  through  the  philosophy  pre- 
viously presented  in  Fosberg  (1972)  where  the 
structural  and  environmental  influences  were  sepa- 
rated. Timelags  obtained  in  the  first  model  were  then 
used  as  a  basis  for  evaluating  environmental  changes 
through  numerical  solutions.  This  model  provided  two 
time-dependent  equations  to  be  used  with  the  similar- 
ity solution  of  the  first  model.  Two  coefficients  of  en- 
vironmental infiuence  were  also  defined  covering  heat 
and  vapor  fiux. 


and 


2x 


[lib] 


Idealized  Analytical  Solutions 


where  the  exponent  x  is  defined  by  the  implicit  relation 


=  1  -  (1  -  (j))^ 


[12] 


If  equation  [12]  is  solved  for  x  by  relaxation,  then  the 
vapor  diffusivity  v  and  thermal  diffusivity  k  in  the  voids 
are  defined  by  the  values  they  would  have  in  the  free 
air.  These  free  air  values  are  indicated  by  the 
subscript . 

These  predictive  and  diagnostic  equations  are  high- 
ly interdependent  and  must  be  solved  simultaneously. 
Because  of  the  nonlinear  character  of  the  basic  void 
diffusion  equations,  they  must  be  solved  numerically 
in  order  to  obtain  valid  prediction  of  duff  and  litter 
moisture. 

The  procedures  required  to  produce  meaningful 
solutions  from  these  equations  require  a  number  of  sep- 
arate solutions  which  are  combined  at  various  stages  to 
provide  the  general  models  and  the  associated  coeffi- 
cients. Two  distinct  models  were  developed.  The  first 
was  for  heat  and  vapor  flux  under  constant  boundary 
conditions  as  might  be  carried  out  in  a  laboratory.  This 
model  was  obtained  first  by  developing  simple  ana- 
lytical solution.  This  analytical  solution  was  then  re- 
fined by  a  series  of  numerical  solutions  to  refine  the 
analytical  coefficients.  Similarity  solutions  were  then 
developed  which  predicted  the  laboratory  timelags  from 
the  physical  properties.  The  numerical  values  used  to 

*Blackmarr,  W.  H.  November  27,  1972,  personal  cor- 
respondence. 


First,  a  simplified  set  of  equations  was  solved 
analytically.  These  analytical  solutions  are  gross 
simplifications  of  the  real  world  situation.  They  are 
useful  in  that  they  provide  a  framework  in  which  to 
evaluate  the  numerical  solutions,  and  they  define  the 
types  of  numerical  simulations  that  are  required  to  ob- 
tain a  prediction  system  of  equations. 

The  analytical  solutions  were  derived  from  sim- 
plified forms  of  equations  [1]  and  [4].  In  addition,  the 
source  and  sink  of  energy  or  mass  were  neglected. 
Thus  the  analytical  solutions  represent  Fickian 
diffusion  through  an  inert  porous  media — one  in  which 
the  particles  respond  instantaneously  to  pore  micro- 
environmental  changes.  With  these  initial  assump- 
tions, the  general  form  of  equations  [1]  and  [4] 
reduced  to 


[13] 


Here,  q  represents  either  the  vapor  density  of  equation 
[1]  or  the  temperature  of  equation  [4],  and  D  is 
either  the  vapor  or  thermal  diffusivity.  If  a  non- 
dimensional  variable  R  is  defined  as 


R(t,z) 


q(t,z)  -  qp.(z) 


then  R  lies  between  0  and  1.  The  variable  q^iz)  repre- 
sents the  final  equilbrium  state  and  q^  is  the  initial 
homogeneous  value  of  q  at  time  0. 
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The  diffusion  equation  [13]  becomes: 


Equation  [17]  then  becomes 


M=  D 
3t 


3^R  _  9R   I  3z 

9z2 


3z    \  (qo-q^) 


9z2 


(R-1) 
(qo-q^) 


[14] 


when  expressed  in  nondimensional  forms  of  q.  The  first 
derivative  of  R  with  respect  to  z  may  be  eliminated  by 
letting 


R  =  Q  exp 


(qo-q„) 


1_  3T 
TD  9t 


1  ,  ,2 

Z  ^ 


-  a 


or 


1 11 
T  9t 


=  -  D 


and 


[18] 


With  this  substitution,  equation  [14]  becomes 


9R 
3t 


gfR 

9z2 


+  R 


9z' 


[1  -  a] 


[15] 


(qo-qe) 


(qo-qe) 


when  a  is  defined  as 


/         iae  \ 

\2  7       (qo-qe)  / 


a  =      (qo-q^)  exp  \± 


If  this  idealized  solution  is  assumed  to  have  fixed 
boundary  conditions  for  q  at  the  top  and  bottom,  then 
the  qg  profile  is  linear  and  the  second  derivatives  of 
qg  vanish. 

This  assumption  then  allows  equation  [15]  to 
be  reduced  further  to 


9R 
9t 


=  D 


9^R 


+  R  T- 


it) 


(qo-q^) 


[16] 


Equation  [16]  may  now  be  solved  by  separation  of 
variables.  We  let  R  equal  the  product  of  two  functions, 
one  dependent  only  on  time  and  one  only  on  depth, 

R  =  T(t)  Z(z) 
so  that  equation  [16]  is 


Z  9T 
D  9t 


where 


T  t4  +  b'TZ 


M 

(qo-qe) 


[17] 


9^Z 


+  (b^  +  a^)  Z  =  0 


[19] 


The  solutions  of  equations  [18]  and  [19]  are 


T  =  exp  (-a^Dt) 


and 


Z  =  A  exp  ((b+ia)z)  +  B  exp  ((b-ia)z) 
so  that 

R  =  TZ  =  exp(-a^Dt) [Aexp((b+ia)z)  +  Bexp (  (b-ia) z) ] 
The  real  part  of  this  solution  is 

R  =  exp(-a^Dt)  exp(bz)  [Ci  cos  az  +  Cz  sin  otz]  [20] 

The  boundary  and  initial  conditions  imposed  on  equa- 
tions [20]  are 

R(0,z)  =  1 
R(t,0)  =  0 
R(t,H)  =  0 

First  applying  the  boundary  condition  at  z  =  0 
0  =  exp(-aDt)C2 


and  at 


thus 


z  =  H 

0  =  exp(-aDt)  exp(bH)  sin  a 
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and 
R  = 


°°  /       2    2  \ 

n=l  ^  ' 

The  general  form  of  the  diffusion  processes  is  obtained 
by  integi-ating  R  over  the  depth  of  the  layer  to  obtain 
the  integral  mean  value 

R=-y  Rdz 


[22] 


OO 

This  integration  shows  that  an  exponential  approach 
to  equilibrium  should  be  expected.  Thus  a  simplified 
form  of  equation  [22]  is 


r;  exp 


[23] 


where  Tg  is  timelag  of  the  diffusion  process.  The  con- 
stant i,  equals  1  since  at  time  zero  R  equals  1.  This 
analytical  solution  applies  only  to  conditions  in  which 
the  external  environment  is  held  constant  and  the  in- 
ternal dynamic  properties  are  initially  uniform.  This 
solution  provides  the  basis  for  evaluation  of  the  system 
dynamics  as  they  would  be  investigated  in  the 
laboratory. 

The  laboratory  type  solution  is  difficult  to  apply  to 
the  field  situation  since  it  requires  constant  boundary 
conditions.  Therefore  a  second  analytical  solution  is 
derived  which  defines  the  functional  form  of  the  dyna- 
mics with  nonstationary  boundary  conditions  and  a 
nonuniform  initial  condition.  This  solution  is  obtained 
by  allowing  the  initial  quantity  to  be  a  function  of  z. 
The  boundary  conditions  are  held  constant  over  some 
small  arbitrary  interval  of  time.  The  value  of  q  from 
the  previous  solution  then  becomes  the  q  for  the  new 
solution,  and  a  new  set  of  boundary  conditions  is 
applied. 

The  solution  for  this  situation  proceeds  along  the 
same  lines  as  that  for  the  homogeneous  case.  The  non- 
dimensional  variable  R  is  now  defined  as: 


The  boundary  and  initial  condition  for  R  from 
solution  of  equation  [20]  are  still  valid.  Thus  the  final 
solution  is  identical  in  form  to  equation  [23]. 

^  /  -St' 

R  =  ^  exp  I  — 


It  is  desirable  to  use  the  rg  defined  in  equation 
[23]  in  this  solution  because  7g  can  be  determined  in 
the  laboratory,  and  the  influences  of  the  physical 
properties  of  the  porous  media  can  be  separated  from 
the  influence  of  the  environment. 
Thus,  if  Tg  is  used: 


R  =  ?  exp 


[24] 


In  this  case,  ^  is  not  necessarily  equal  to  1  during 
the  moisture  exchange  since  it  expresses  the  degree  of 
nonhomogeneity  when  a  new  set  of  boundary  con- 
ditions is  applied. 

The  forms  of  equations  [23]  and  [24]  provide 
the  basis  for  evaluating  the  numerical  simulation. 
Neither  equation  provides  the  necessary  detail  to  pre- 
dict theoretically  the  operational  coefficients  required 
for  vapor  density  or  temperature  flux. 

The  analytical  solution  for  fixed  boundary  con- 
dition [equation  23]  suggests  that  numerical  solu- 
tions and  laboratory  experiments  involving  struc- 
turally homogeneous  layers  of  litter  and  duff  are  re- 
quired to  define  the  timelags.  The  equations  suggest 
that  an  exponential  approach  to  equilibrium  is  to  be 
expected.  Since  the  initial  vapor  or  temperature  pro- 
files are  not  homogeneous,  a  different  response  time 
should  be  expected  for  each  different  set  of  initial 
conditions.  Thus,  if  the  timelags  as  derived  in  the  lab- 
oratory are  to  be  extended  to  field  applications,  then  a 
second  coefficient  expressing  this  nonhomogeneous 
initial  state  would  be  required.  This  step  can  be  ac- 
complished by  either  allowing  the  timelag  to  be 
environmentally  dependent,  or  by  allowing  the  new 
coefficient,  ^,  to  be  variable.  This  paper  is  based  on 
the  more  general  solutions  of  a  variable  homogeneity 
coefficient  with  laboratory-derived  timelags. 


(qo(z)  -  qe(z)) 


and  is  substituted  into  equation  [13].  As  in  the  pre- 
vious solution,  R  is  redefined  as 


2y  (qo-qe) 


to  eliminate  the  first  derivative.  Again  the  general 
solution  is  obtained  by  separation  of  variables.  In  this 
case,  the  integral  for  time  is  taken  over  8t,  an  arbitrary 
time  interval. 


Relation  of  Timelag  to  Physical  Properties 


The  analytical  solutions  presented  above  predicted 
an  exponential  approach  to  the  equilibrium  profile.  If 
this  process  is  valid  in  the  more  complex  physical 
system  found  in  nature,  then  the  gain  or  loss  of  tem- 
perature and  vapor  density  can  be  expressed  through 
timelag  determination.  This  timelag  is  directly  related 
to  the  physical  structure  of  litter  and  duff  layers.  The 
analytical  solutions  showed  that  the  coefficients  in  the 
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exponent  were  related  to  the  physical  variables  of  layer 
thickness  and  bulk  diffusivity. 

If  nondimensional  groups  of  timelag,  bulk  or 
simple  diffusivity,  and  thickness  are  formed,  then  a 
similju-ity  number  is  defined  as: 


may  be  similarly  expressed  tor  heat  and  vapor  transfer. 
The  Fourier  number  for  heat  transfer  is: 

Ft  =  ^  [28] 


^1  "  H2 


and  for  water  vapor 


[25] 


which  is  readily  recognized  as  the  Fourier  number 
(Gukhman  1965). 

The  complete  model  proposed  in  this  paper  in- 
cludes a  number  of  other  variables.  These  are  the 
diffusivity  in  the  voids,  and  the  particle  timelag  as  well 
as  the  depth  and  horizon  timelags.  Also,  since  the  part- 
icles exchange  heat  and  vapor  with  the  voids,  the  total 
mass  of  particles  and  their  volume  should  also  be 
considered.  Three  nondimensional  groups  of  variables 
may  be  defined.  The  ratio  of  the  bulk  density  to  part- 
icle density — the  packing  ratio  — is  one  minus  the 
porosity.  The  second  group  is  the  ratio  of  the  system 
timelag  to  particle  timelag.  The  third  is  composed 
of  thickness  of  the  system,  the  diffusivities  in  the  voids, 
and  a  time  variable.  This  time  variable  was  chosen  as 
the  particle  timelag  since  the  descriptive  model 
suggested  this  time  interval  as  the  critical  time 
variable.  The  three  groups  are: 


(6) 


for  mass  exchange  and 


-fe)'(^) 


for  heat  exchange.  These  variables  form  similarity 
numbers  of: 


and 


F„  = 


fe)(^) 


[26] 


[27] 


The  packing  ratio  has  been  placed  in  the  denom- 
inator of  equations  [26]  and  [27],  since  the  timelag 
should  increase  as  the  bed  porosity  decreases.  The 
Fourier  number  of  the  analytical  solutions  [equation  25] 


BM  V 
H2 


[29] 


The  hat  (  )  over  the  diffusivities  in  equations  [28]  and 
[29]  denote  the  net  value  of  the  diffusivity  obtained 
by  neglecting  the  particle  effect  in  diffusion  experi- 
ments. 

If  the  detailed  theoretical  results  are  a  more 
general  solution  of  the  physical  processes  than  those 
from  simple  diffusions,  then  the  timelags  predicted 
by  either  method  should  be  the  same  for  a  given  set  of 
physical  conditions. 


Numerical  Solution  Methods 

The  analytical  and  similarity  functions  provided 
the  basis  for  evaluating  numerical  solutions  and 
observational  data.  The  complete  set  of  equations 
[1  through  6]  must  be  solved  simultaneously.  Equa- 
tions 1  and  4  are  highly  nonlinear  and  therefore  must 
be  solved  numerically  if  the  full  dynamics  are  to  be 
preserved. 

Finite  difference  schemes  must  be  adopted  for  both 
the  spatial  and  temperal  derivatives  to  solve  these 
equations.  The  philosophy  in  describing  the  space 
derivatives  was  that  certain  regions  of  the  system 
would  require  very  high  resolutions,  while  in  others, 
the  dynamic  process  would  be  slow  and  would  not  be 
characterized  by  strong  gradients.  Techniques  utiliz- 
ing a  constant  space  interval  would  require  the  same 
high  resolution  in  the  quiescent  regions  as  that  desired 
in  the  strong-gradient,  rapid-change  regions.  The 
finite  difference  scheme  adopted  allowed  for  arbi- 
trarily selected  space  intervals  between  each  com- 
putational point.  This  preserved  the  high  resolution 
desired  in  certain  regions,  and  reduced  the  total  effort 
required  to  obtain  a  particular  solution.  The  particular 
finite  difference  stencils  adopted  for  the  void  to  void 
diffusion  are  based  on  the  mean  value  between  each 
computational  point.  Thus,  the  first  derivative  of  a 
scalar  quantity  (S)  at  a  point  i  is: 


S.^,  -H  S.      S.  -(-  S.  , 

1+1  1  X  1-1 


Z.  ^  +  Z.  Z.  +  Z.^, 
1-1        i        1  1-1-1 
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where  Z  is  the  depth  below  the  surface.  If  the  computa- 
tional points  are  equally  spaced,  then  this  scheme  re- 
duces to  the  familiar  centered  difference  method.  If  the 
computational  points  are  not  equally  spaced  and  the 
redundant  value  at  i  is  eliminated,  then  the  derivative 
is  not  defined  at  the  point  i.  The  second  derivative  is 
defined  by  the  difference  between  the  uncentered  first 
derivatives 

and 

so  that  the  second  derivative  is: 

2  ^  / 

This  derivative  also  reduces  to  the  familiar  centered 
second  derivative  with  equally  spaced  computational 
points. 

The  finite  difference  time  step  used  in  this  analysis 
was  the  Adams-Bashforth  Scheme  (Lilly  1965): 


where  the  subscripts  t+1,  t,  and  t-1  refer  to  the  new 
time,  current  time,  and  previous  time,  separated  by 
the  time  interval  8t.  This  scheme  was  chosen  because 
of  well  defined  stability  in  nonlinear  thermofluid 
dynamic  solutions. 

The  time  step  8t  must  be  made  sufficiently  small 
so  that  the  characters  of  the  continuous  equations  are 
preserved  with  a  minimum  of  distortion.  The  maximum 
value  of  8t  is  then  defined  by  linear  stability  criteria  of 

max      2D . /  . 

i  min 

where  8z  is  the  distance  between  computational  points 
and  D  is  the  diffusivity  of  heat  or  vapor  in  the  void  to 
void  exchange. 

A  second  constraint  not  related  to  computational 
stability,  must  be  placed  on  the  magnitude  of  the  time 
step.  This  constraint  arises  because  of  the  large 


differences  between  the  internal  diffusivities  of  the 
particles  and  those  in  the  voids.  This  constraint  is  not 
a  mathematical  or  physical  constraint,  but  is  directly 
related  to  the  resolution  obtainable  in  digital  com- 
puters and  to  the  range  over  which  the  particle  homo- 
geneity coefficient  is  properly  defined.  If  the  two  dif- 
ferential equations  describing  the  heat  and  mass 
processes  in  the  particles  are  included  with  those 
describing  bed  processes,  then  the  number  of  compu- 
tations that  must  be  made  becomes  extremely  large 
and  beyond  capabilities  of  even  the  largest  computers. 
Computational  limits  are  imposed  by  the  modified 
analytical  solution  [equations  3  and  6].  This  restriction 
is  covered  by  i,  defined  only  over  a  small  range  of 
St 

The  net  effect  is  to  increase  stability.  This  tradeoff 
increases  the  number  of  time  steps,  but  drastically 
decreases  the  total  number  of  computations  to  the  point 
that  a  net  gain  is  obtained  and  practical  solutions 
are  possible. 

Results  of  the  Computations 

The  physical  properties  of  litter  and  duff  required 
to  carry  out  the  numerical  simulations  were  the  depth 
of  the  horizon,  bulk  density,  particle  density,  thermal 
and  moisture  timelags  of  the  particles,  specific  heat 
of  the  particles,  and  the  equilibrium  moisture 
isotherms. 

The  two  species  for  which  a  complete  analysis  was 
carried  out  were  ponderosa  pine  (Pinus  ponderosa)  and 
lodgepole  pine  (Pinus  contorta).  The  particle  densities 
for  ponderosa  pine  were  obtained  from  Brown  (1970), 
and  for  lodgepole  pine  through  personal  correspondence 
with  him.^  The  moisture  timelags  and  equilibrium 
moisture  isotherms  of  both  species  for  freshly  fallen  and 
weathered  needles  were  obtained  from  Blackmarr 
through  personal  correspondence  or  unpublished  re- 
search (see  footnote  4).  The  thermal  timelags  were 
calculated  from  the  relations  developed  previously 
(Fosberg  1973).  The  specific  heat  was  given  by  equation 
[10].  The  bulk  densities  of  the  horizons  were  obtained 
from  samples  collected  on  the  Manitou  and  Fraser 
Experimental  Forests  of  the  Rocky  Mountain  Forest 
and  Range  Experiment  Station,  and  on  the  Pingree 
Park  Campus  of  Colorado  State  University.  The  pon- 
derosa pine  data  came  only  from  the  Manitou  Experi- 
mental Forest,  while  the  lodgepole  pine  data  came  from 
all  three  areas.  Data  used  in  the  simulations  are 
summarized  in  table  1. 

These  data  are  supplemented  through  arbitrarily 
chosen  variations  for  each  physical  property.  Each  of 
three  physical  variables — depths,  bulk  density,  and 
particle  timelag — were  systematically  varied  over  a 
large  range  of  conditions  while  the  others  were  held 
constant.  The  depths  were  selected  as  1,  3,  4,  6,  and  8 
cm.  The  selected  particle  timelags  were  3,600 
seconds,  7,200  seconds,  and  the  measured  4,608 
seconds.  While  bulk  densities  of  litter  are  mostly 

^Brown,  J.  K.  October  11, 1972 ,  personal  correspondence. 
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Table  1. --Physical  properties  of  duff  and  litter  horizons 


Horizon 


Particle  timelags  (sec.) 


Density  (gm  cm"^) 


Mo 


isture  (t^)     Thermal   (t^)        Particle  (pp)     Bed  bulk  (pg) 


Li  tter  Aoo 

Fresh  1 y  fa  1 len 
needles 

Weathered 
need  1 es 

Duff  Ao 


Litter  Ac 

Freshly  fal len 
needles 

Weathered 
needles 


Duffo 


l'<,760 


i»,608 


63,000 


3,636 


Ponderosa  Pine 


0.  12 


11 


0.51 
.h7 


Lodgepole  Pine 


.10 


10 


.56 

.56 
.56 


0.08 

.08 
.26 


.09 

.09 
.36 


in  the  range  of  0.06  to  0.1  gm/cm^,  the  mean  value 
of  0.08  gm/cm^  was  supplemented  with  bulk  den- 
sities 0.01,  0.02,  and  0.03  gm/cm^  to  provide  a 
larger  range  of  values  to  evaluate  the  similarity 
numbers. 

The  individual  computer  simulations  using 
equations  [1]  through  [12]  for  a  particular  value  of 
particle  and  bed  properties  were  evaluated  through  the 
characteristics  provided  by  the  analytical  solutions 
[equations  23  and  24],  and  by  the  similarity  relations 
predicted  by  equations  [26]  and  [29]. 

The  first  set  of  simulation  results  applies  to  a 
constant  environment.  These  simulated  results  de- 
fined the  timelag  characteristics  and  provided  the 
functions  to  predict  them  from  observations  of  the 
physical  structure. 

The  second  set  of  simulation  results  utilizes  the 
predicted  constant  responses  and  applies  them  to 
prediction  of  heat  and  moisture  response  with  time- 
dependent  boundary  conditions.  This  second  set  of 
results  provides  the  framework  for  field  prediction. 

Constant  Boundary  Conditions 

The  numerical  simulations  of  these  conditions 
were  conducted  under  the  conditions  of  an  initially 


homogeneous  moisture  and  temperature.  At  an  arbi- 
trary time  zero,  an  arbitrary  step  change  of  either  rel- 
ative humidity  or  temperature  was  imposed  while  the 
other  boundary  was  held  at  its  original  value.  The 
step  change  in  the  boundary  condition  then  set  up  a 
time-dependent  process  from  the  initial  homogeneous 
state  to  one  in  which  the  final  equilibrium  profile  was  a 
linear  function  with  depth.  The  time-dependent 
approach  to  the  equilibrium  could  then  be  used  to 
determine  the  timelag  of  the  complete  system.  As 
examples,  several  of  the  simulations  for  moisture 
exchange  are  shown  here  to  illustrate  the  details  of 
the  nonlinear  solutions. 

The  first  set  of  examples  is  for  isothermal  con- 
ditions. Boundary  conditions  were  set  at  20  percent 
relative  humidity  at  the  surface,  and  the  initial  homo- 
geneous void  relative  humidity  was  80  percent.  The 
temperature  was  300  K.  Particle  and  bed  properties 
were  defined  by  weathered  ponderosa  pine  litter.  The 
depths  of  the  two  horizons  were  arbitrarily  selected  as 
3  and  4  cm  (reasonable  values  for  stands  of  basal  area 
in  the  15  to  20  m^/hectare  range  along  Colorado's  east 
slope). 

Initial  moisture  properties  were  used  to  evaluate 
the  timelags.  The  integral  particle  moisture  content  for 
the  two  examples  (fig.  2)  follows  the  expected  expo- 
nential decay.  The  3-cm-thick  layer  had  a  timelag  of 
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Figure  2. — Time  response  of  integral  particle  mois- 
ture content  in  a  homogeneous  layer  of  weathered 
ponderosa  pine  litter  under  constant  boundary 
conditions.  Bulk  density  is  0.08  gm  cm  I  Particle 
timelag  is  4,608  sec.  Curve  A  is  for  a  layer  3  cm 
thick.  Curve  B  is  for  a  layer  4  cm  thick. 
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4,248  seconds  while  the  4-cm  layer  had  a  timelag  of 
7,560  seconds.  The  integral  of  void  vapor  density 
(fig.  3)  for  the  two  cases  also  showed  the  exponential 


decay,  and  had  timelags  nearly  identical  to  those  for 
particle  moisture  content.  Void  vapor  density  was 
chosen  to  represent  the  moisture  transport  since  the 
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Figure  3. — Time  response  of  integral  void  vapor  den- 
sity in  a  homogeneous  layer  of  weathered  pon- 
derosa pine  litter  under  constant  boundary 
conditions. Bulk  density  is  0.08  gm  cm~^.  Particle 
timelag  is  4,608  sec.  Curve  A  is  for  a  layer  3  cm 
thick.  Curve  B  is  for  a  layer  4  cm  thick. 
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equilibrium  moisture  content  of  the  particles  is  a  non- 
linear function  of  relative  humidity,  and  because 
relative  humidity  is  strongly  dependent  on  tempera- 
ture. Void  relative  humidity  profiles  for  the  two 
examples  (fig.  4)  decayed  rapidly  to  the  linear  equi- 
librium profile.  The  particle  moisture  content  profiles 
(fig.  5 )  can  be  predicted  from  the  void  vapor  density  and 


the  temperature,  since  the  vapor  density  and  the  mois- 
ture content  approach  equilibrium  at  the  same  rate. 

Temperature  in  the  duff  and  litter  layer  changed 
much  more  rapidly  than  moisture.  The  integral  temp- 
erature change  over  time  also  approached  equilibrium 
exponentially  (fig.  6).  Surface  boundary  conditions  for 
the  temperature  experiments  were  290,  296,  305,  and 
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Figure  4. — Profiles  of  void  relative  humidity  for  wea- 
thered ponderosa  pine  litter  under  constant 
boundary  conditions.  Bulk  density  is  0.08  gm 
cnn"^  Particle  timelag  is  4,608  sec.  Numbers  on 
lines  indicate  hours  since  beginning  of  solution.  A 
is  for  a  layer  3  cm  thick;  B  is  for  a  layer  4  cm 
thick. 


Figure  5. — Profiles  of  particle  moisture  content  for 
weathered  ponderosa  pine  litter  under  constant 
boundary  conditions.  Bulk  density  is  0.08  gm 
cm"l  Particle  timelag  is  4,608  sec.  Numbers  on 
lines  indicate  hours  since  beginning  of  solution. 
A  is  for  a  layer  3  cm  thick;  B  is  for  a  layer  4  cm 
thick. 
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Figure  6. — Time  response  of  integral  temperature  of 
particles  for  a  homogeneous  litter  layer  under 
constant  boundary  conditions.  Bulk  density  is 
0.03  gm  cm  Layer  is  8  cm  thick.  Particle  time- 
lag  is  6  sec. 


Time (hours) 


310  K.  The  initial  temperatures  were  300  K.  Initial 
and  boundary  condition  relative  humidities  were  20 
percent.  Profiles  of  temperature  (fig.  7)  behaved 
similarly  to  those  for  vapor  density  in  that  their 
equilibrium  profiles  were  linear.  When  both  moisture 
and  temperature  boundary  condition  changes  were 
imposed,  the  vapor  density  profiles  still  decayed  to 
the  linear  equilibrium  profile  (fig.  8).  The  relative 
humidity  boundary  condition  was  20  percent  and  the 
temperature  boundary  condition  was  310  K.  Initial 
conditions  were  80  percent  relative  humidity  and 
300  K.  The  internal  maximum  of  vapor  density  repre- 
sents the  thermally  driven  moisture  flux.  Thus,  the 
conservative  nature  of  temperature  and  vapor  density 
provides  a  more  stable  basis  for  predicting  moisture 
content  than  methods  based  on  relative  humidity. 

The  individual  results  of  the  25  simulations  where 
each  bed  and  particle  property  was  systematically 
varied  were  combined  by  the  similarity  results  defined 
by  equations  [26]  and  [27].  These  similarity  relations 
provided  the  basis  for  predicting  the  horizon  timelags 
from  the  physical  properties  of  particles  and  bed.  The 
vapor  density  relationship  was; 
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and  is  graphically  shown  in  figure  9.  Similarly  the 
relationship  for  heat  transfer  is 
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Figure  7. — Profiles  of  particle  temperature  for  homo- 
geneous litter  layer  under  constant  boundary  con- 
ditions. Bulk  density  is  0.03  gm  cm  l  Layer  is  8 
cm  thick;  particle  timelag  is  6  sec.  Numbers  on 
lines  indicate  hours  since  beginning  of  solution. 
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Figure  8. — Profiles  of  vapor  density  for  homogeneous 
litter  layer  under  constant  boundary  conditions 
with  temperature  and  humidity  stress.  Bulk 
density  is  0.03  gm  cm  I  Particle  moisture  time- 
lag  is  4,608  sec,  particle  temperature  timelag  is  6 
sec.  Layer  is  8  cm  thick.  Numbers  on  lines  indicate 
hours  since  beginning  of  solution. 
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■Response  timelags  of  litter  and  duff  layers. 
Line  1  is  for  bulk  density  of  0.08  gm  cm"^  particle 
timelag  of  4,608  sec.  Line  2  is  for  bulk  density  of 
0.03  gm  cm"^  particle  timelag  of  7,200  sec.  Line  3 
is  for  bulk  density  of  0.02  gm  cm  ^  particle  time- 
lag  of  7,200  sec.  Line  4  is  for  bulk  density  of  0.01 
gm  cm  ^  particle  timelag  of  3,600  sec.  Line  5  is  for 
bulk  density  of  0.01  gm  cm  ^  particle  timelag  of 
7,200  sec. 
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The  similarity  numbers  of  equations  [30]  and  [31] 
are  averages  of  the  individual  simulations.  The 
range  of  the  similarity  numbers  was  20  percent  of  the 
mean  for  mass  transfer  and  10  percent  for  heat  trans- 
fer. These  two  similarity  relations  provided  the  basis 
for  calculating  the  moisture  and  temperature  time- 
lags  under  constant  or  laboratory  conditions. 

Time-Dependent  Boundary  Conditions 

The  natural  environment  rarely  provides  constant 
boundary  conditions.  IMore  typically,  the  boundary 
conditions  are  continuously  changing  temperatures 
and  humidities.  Functions  must  therefore  be  establish- 
ed for  these  variations  under  which  the  horizons  do  not 
reach  an  equilibrium  state.  As  with  the  constant 
environment  case,  the  analytical  solutions  provide  the 
functional  framework  for  evaluation  of  numerical 
solutions.  Laboratory  and  field  solutions  differ  only  by 
a  linear  multiplier  of  the  exponential  terms.  This 
linear  coefficient  reflects  the  inhomogeneity  of 
the  initial  conditions.  The  coefficient  C  of  equation 
[24]  would  be  unity  under  constant  boundary  con- 
ditions, but  will  be  different  from  unity  under  nearly 
all  field  conditions. 

The  time-dependent  boundary  conditions  were 
defined  by  step  changes  of  arbitrary  duration  and 
arbitrary  magnitude.  These  changes  established  a 
nonhomogeneous  interval  profile.  The  integral  vapor 
density  (fig.  10)  and  the  integral  temperature  (fig.  11) 
change  over  time  show  two  marked  characteristics. 
The  responses  of  the  layers  are  much  more  rapid  when 
they  start  from  the  nonhomogeneous  profile  than  when 
they  start  from  the  homogeneous  profile.  This  dif- 
ference can  be  interpreted  either  as  a  change  in  the 
timelag  or  as  enhanced  transport  due  to  larger 
gradients.  Because  the  timelag  was  defined  in  terms  of 
physical  properties,  the  second  interpretation  is  pre- 
ferred. This  latter  interpretation  has  the  advantage 
that  dynamic  and  physical  characteristics  are  sepa- 
rated in  the  functions. 

The  coefficient  C,  of  equation  [24]  then  must  be 
defined  for  the  heat  and  mass  transfers.  This  was 
accomplished  by  calculating  the  actual  change  in  vapor 
density  or  temperature,  the  potential  change  in  vapor 
density  or  temperature,  and  using  the  timelags  defined 
from  the  constant  boundary  condition  simulation.  This 
C  is: 


The  homogeneity  coefficient  for  mass  transfer  is 
0.42  and  for  heat  transfer  is  0.755.  No  variation  of  the 
coeffiicient  was  found  in  the  simulation  used  to 
determine  these  values.  Variation  probably  does  exist 
and  would  be  found  if  more  simulation  or  field  experi- 
ments are  conducted. 


Comparison  of  Generalized  Solutions 
with  Fickian  Solutions 

Bulk  diffusivities  of  the  horizons  can  be  used 
to  validate  the  theoretical  predictions.  These  bulk 
diffusivities  can  be  calculated  from  the  two  pairs  of  sim- 
ilarity equations.  By  solving  equations  [29]  and  [30] 
for  Tg]y[  and  equations  [27]  and  [31]  for  rg-p,  the  bulk 
vapor  diffusivity  v  and  thermal  diffusivity  k  can  be 
calculated.  The  effective  vapor  diffusivity  is: 

_2 

0.18  (H  v„)3 
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and  the  effective  thermal  diffusivity  is: 
2 

0.14  (H  k_,)3 


The  coefficients  0.18  and  0.14  are  taken  from  Fosberg 
(1972,  1973). 

The  vapor  diffusivity  for  the  example  in  figures  2 
through  5  is: 


3.1  X  10  ^  cm^/sec. 


The  thermal  diffusivity  for  the  temperature  change  of 
figure  6  is: 


5.7  X  10  3  cm^/sec. 


exp  -  — 

where  6S  is  the  actual  change  of  the  scalar  quantity, 
AS  is  the  potential  change,  and  t^  is  the  duration  of  the 
pulse. 


These  bulk  diffusivities  compare  well  with  experiment- 
al values  based  on  equilibrium  solutions  (Lettau  1954; 
Jackson  1964a,  19646,  1964c  ). 

Reported  moisture  response  rates  for  litter  and 
duff  have  typically  been  very  slow.  Johnson's  (1968) 
analysis  of  field  data  gives  timelags  an  order  of  magni- 
tude greater  than  reported  here.  Since  Van  Wagner 
(1970)  used  the  timelag  as  the  environment  dependent 
variable,  it  is  difficult  to  determine  the  exact  magni- 
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Figure  10. — Time  response  of  integral  void  vapor  den- 
sity for  homogeneous  litter  layers  under  varying 
boundary  conditions.  Bulk  density  is  0.03  gm  cm  ^ 
Particle  timelag  is  3,600  sec.  Thickness  is  8  cm.  Lines 
A  and  B  had  a  20  percent  relative  humidity  on  the 
top  surface  and  a  40  percent  relative  humidity  on 
the  bottom  surface.  The  top  surface  was  changed  to 
80  percent  relative  humidity.  Line  C  had  a  20 
percent  relative  humidity  on  the  top  surface  and 
an  80  percent  relative  humidity  on  the  bottom 
surface.  The  top  surface  was  changed  to  40  percent 
relative  humidity. 


Time  ( hours 


tude  of  the  structure  timelag.  Examination  of  the  forc- 
ing and  response  functions  suggests  a  magnitude 
similar  to  that  reported  by  Johnson  (1968). 

Data  used  by  Johnson  (1968)  and  Van  Wagner 
(1970)  come  from  observations  taken  once  a  day.  Thus, 
the  reported  responses  were  artificially  forced  to  long 
timelags  since  the  minimum  resolvable  timelag  cannot 
be  less  than  1  day.  Data  collected  at  high  sample  rates^ 
show  high  response  rates  and  short  timelags  similar 
to  those  predicted  by  the  theoretical  solutions.  Thus, 
discrepancies  between  the  theoretical  solution  and 
timelag  values  previously  reported  are  due  to  artificial 
filtering  created  through  sampling  procedures. 

Comparisons  of  the  theoretical  prediction  to 
observation  were  made  with  Nelson's  (1969)  experi- 
mental results  in  sawdust  layers.  The  theory  predicted 
a  timelag  of  360  minutes,  while  Nelson  observed  time- 
lags  that  ranged  from  390  to  420  minutes  for  low 
moisture  timelags.  Measured  timelags  of  ponderosa 

^McLeod,  Bruce.  December  7,  1973,  Microwaves,  a  new 
tool  for  forest  and  watershed  management.  15  p.  (Unpublished 
report,  Mont.  State  Univ.,  Bozeman.J 
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Figure  11. — Time  response  of  integral  particle  temper- 
ature for  homogeneous  litter  layer  under  varying 
boundary  conditions.  Bulk  density  is  0.03  gm 
cm  I  Particle  timelag  is  6  sec.  Layer  is  8  cm  thick. 
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pine  needles  in  loosely  packed  beds''  were  close  to  those 
predicted  by  the  theoretical  solutions. 

Conclusions 

The  theoretical  model  developed  in  this  paper  is  a 
generalization  of  previously  derived  theoretical  solu- 
tions. The  improvement  over  traditional  solutions  was 
obtained  by  following  the  physical  processes  more 
closely.  Three  main  quantitative  results  were  obtained 
from  the  theoretical  predictions. 

The  first  was  the  prediction  of  timelags  from 
physical  properties  of  particles  and  from  horizon 
structure.  These  thermal  and  vapor  timelags  provide 
basic  coefficients  for  analytical  predictions  of  heat  and 
vapor  exchanges.  Thermal  and  vapor  timelags  cannot 
be  calculated  as  the  particle  timelags  go  to  zero.  Non- 
sorbing  porous  media  flow  is  depicted  by  the  simple  ana- 
lytical solution  and  similarity  is  not  preserved.  While 
this  limitation  indicates  that  the  solutions  are  not 
completely  generalized,  calculations  over  the  range  of 
properties  in  organic  soils  are  valid. 

The  second  primary  result  was  two  pairs  of  pre- 
diction equations.  The  first  pair  apply  to  heat  and  vapor 
changes  with  constant  boundary  conditions  and  ini- 
tially homogeneous  profiles.  These  equations  depict 
the  behavior  of  the  horizons  in  the  laboratory.  The 
second  pair  of  equations  for  heat  and  vapor  change  pre- 
diction relax  the  requirement  of  initially  homogeneous 
profiles  and  constant  boundary  conditions.  This  pair  of 
equations  depicts  the  processes  under  natural  or  field 
conditions  of  temperature  and  moisture  profiles. 

These  solutions  for  field  conditions  permit  a  vari- 
able set  of  boundary  conditions.  The  equation  requires 
surface  temperature  and  the  vapor  density  at  the  sur- 
face. Because  the  boundary  conditions  are  held 
constant  over  a  small  integration  interval,  but  can  vary 
between  finite  intervals,  the  boundary  conditions  are 
actually  finite  difference  approximations  to  the  con- 
tinuous change.  The  homogeneity  coefficient  in  these 
field-applicable  equations  reflects  the  departure  of  the 
initial  state  from  homogeneous  temperature  and 
moisture  profiles.  If  most  of  the  changes  in  temperature 
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and  moisture  are  determined  by  surface  boundary 
conditions,  then  the  homogeneity  coefficient  will  be 
less  than  one.  A  coefficient  of  less  than  one  then  indi- 
cates a  deep,  slowly  changing  reservoir.  A  coefficient 
of  one  applies  only  in  the  laboratory.  Coefficients 
greater  than  one  indicate  quiescent  surface  conditions, 
with  most  of  the  change  taking  place  in  the  deeper 
horizons. 

The  third  quantitative  result  was  the  recalculation 
of  the  bulk  diffusivities  for  nonsorbing  particles.  These 
bulk  diffusivities  are  the  same  as  those  traditionally 
determined  in  the  laboratory,  and  are  based  on  the 
assumption  of  instantaneous  response  of  the  particles. 

The  model  requires  horizon  thickness,  porosities, 
and  particle  timelags  to  make  the  predictions.  The 
particle  timelags  must  be  determined  in  the  laboratory. 
Neither  particle  timelags  nor  densities  vary  much 
between  samples.  Bulk  densities  of  each  horizon  also 
vary  within  a  narrow  range,  so  the  porosity  can  be 
calculated  from  the  particle  and  bulk  densities. 
Horizon  depth  varies  widely  even  in  small  plots  of 
constant  basal  area,  so  that  there  is  a  large  potential 
uncertainty  in  extrapolation  of  an  individual  calcu- 
lation to  different  sites.  The  only  other  data  require- 
ments of  the  model  are  the  boundary  conditions.  These 
data  are  the  temperature  and  atmospheric  moisture 
content  at  the  upper  surface. 

The  model  can  be  applied  to  forestry  problems  such 
as  fire  danger  by  specifying  the  physical  properties  of 
thickness,  particle  timelag,  and  porosity.  With  these 
data  the  vapor  and  thermal  timelags  can  be  calculated. 
The  field  prediction  equations  are  then  solved  sequen- 
tially for  each  horizon.  The  external  boundary  condi- 
tions of  surface  temperature  and  vapor  density  apply  to 
the  surface  horizon.  The  temperature  and  vapor  den- 
sity boundary  conditions  for  internal  horizons  are 
specified  by  the  vapor  density  and  temperature  of  the 
adjacent  higher  horizon.  The  lowest  horizon  must  be 
selected  sufficiently  deep — such  as  below  the  depth  of 
diurnal  fluctuation  for  prediction  extending  a  few  days 
or  below  the  depth  of  annual  fluctuation  for  seasonal 
calculation — to  eliminate  the  lower  boundary  con- 
dition. Since  the  moisture  content  timelag  and  vapor 
density  timelag  are  nearly  identical,  the  moisture  con- 
tent profiles  can  be  calculated  directly  from  the 
equilibrium  moisture  isotherms.  These  theoretical 
solutions  do  not  contain  liquid  water  processes.  Addi- 
tion of  liquid  water  can  be  considered  through  minor 
modification  of  classic  infiltration  theory  and  therefore 
complete  the  model  of  practical  application. 
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Appendix — List  of  Symbols 


Symbol  Definition  Units 

Cp  Specific  heat  of  air  constant  pressure  cal  gm  ^  °K  ^ 

.2 

e  Vapor  pressure  in  air  dynes  cm 

.2 

Bg  Saturation  vapor  pressure  dynes  cm 

Similarity  number  for  moisture  transfer  in  litter 

and  duff  dimensionless 

f'j^^  Similarity  number  for  moisture  transfer  in 

nonsorbing  porous  media  dimensionless 

Frjy  Similarity  number  for  heat  transfer  in  litter  and 

duff  dimensionless 

t'l  Similarity  number  for  heat  transfer  in 

nonsorbing  porous  media  dimensionless 

H  Thickness  of  a  homogeneous  horizon  cm 

h  Relative  humidity  percent 

2  -1 

kgjk  Thermal  diffusivity  in  voids  cm  sec 

2 

kg  Thermal  diffusivity  in  free  air  cm  sec 

2  -1 

k  Effective  thermal  diffusivity  in  litter  or  duff  cm  sec 

Mass  of  air  in  litter  or  duff  volume  gm 

Mp  Mass  of  litter  particles  in  litter  or  duff  volume  gm 

m  Moisture  content  of  litter  or  duff  (mass  of 

contained  water  to  ovendry  mass  of  particle)  percent 

mg  Equilibrium  moisture  content  percent 

mj  Initial  moisture  content  percent 

q  Arbitrary  scalar  variable  dimensionless 

qg  Equilibrium  value  of  arbitrary  scalar  dimensionless 

q^j  Initial  value  of  arbitrary  scalar  dimensionless 

R  Gas  constant  for  water  vapor  cal  gm  ^  °K  ^ 

s  Source  or  sink  of  heat  in  litter  or  duff  cal 

T  Temperature  in  litter  or  duff  deg  K 

t  Time  coordinate  sec 
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Appendix — List  of  Symbols  (Continued) 


Symbol  Definition  Units 

Vg  Volume  of  litter  or  duff  cm 

Vp  Volume  of  particles  in  litter  or  duff  cm 

Vy  Volume  of  voids  in  litter  or  duff  cm 

z  Vertical  space  coordinate  cm 

/3  Packing  ratio  (1  -  porosity)  dimensionless 

Am  Potential  change  of  particle  moisture  content  percent 

AT  Potential  change  of  particle  temperature  deg  K 

8H^  Change  of  heat  in  litter  or  duff  voids  cal 

6Hp  Change  of  heat  in  litter  or  duff  particles  cal 

5m  Change  of  moisture  content  in  particles  percent 

6Tp,8T  Change  of  temperature  in  particles  deg  K 

8t  Finite  time  step  sec 

Homogeneity  coefficient  for  moisture  transfer  dimensionless 

Homogeneity  coefficient  for  heat  transfer  dimensionless 

.3 

Pa  Density  of  air  gm  cm 

.3 

Pb  Bulk  density  of  litter  or  duff  gm  cm 

,3 

Pf  Density  of  litter  or  duff  particles  gm  cm 

.3 

Pv  Water  vapor  density  in  voids  gm  cm 

Pv'  Water  vapor  density  equivalent  of  moisture  in  ^ 

particles  gm  cm 

Tgjyj  Moisture  timelag  of  homogeneous  horizon  of 

litter  or  duff  sec 

Tg-p  Thermal  timelag  of  homogeneous  horizon  of  litter 

or  duff  sec 

Tpjyj  Moisture  timelag  of  litter  or  duff  particle  sec 

Tpp  Thermal  timelag  of  litter  or  duff  particle  sec 

V,  I'g  Moisture  diffusivity  in  voids  cm^  sec  ^ 

Moisture  diffusivity  in  free  air  cm^  sec  ^ 
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Appendix — List  of  Symbols  (Continued) 


Symbol  Definition  Units 

2  -1 

V  Effective  moisture  diffusivity  in  litter  or  duff  cm  sec 

(b  Porosity  of  litter  or  duff  (Volume  of  voids/volume 

of  horizon)  dimensionless 

X  Source  or  sink  of  moisture  in  litter  or  duff  percent 
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